% This function estimates the score of interest, which matters for later
% covariance matrix estimation.
function [S] = score_estimation(method,data)

switch lower(method)
    case {'mean','ols'}
        
        D = data.D;
        X = data.X;
        Y = data.Y;
        
        MD = D-X*(X\D);
        X_mat = [D X];
        beta_hat = X_mat\Y;
        U_hat = Y-X_mat*beta_hat;
        S = MD.*U_hat/mean(MD.^2); % score
        
    case 'iv'
        
        D = data.D;
        X = data.X;
        Y = data.Y;
        Z = data.Z;
        n = length(D);
        
        X_mat = [D X];
        Z_mat = [Z X];
        beta_iv = (Z_mat'*X_mat)\(Z_mat'*Y);
        U_hat = Y-X_mat*beta_iv;

        % In IV cases, beta0_iv-beta0=(Z1'*M2*X1)^{-1}(Z1*M2*U)
        M2 = eye(n)-X*((X'*X)\X');
        S = (Z'*M2)'.*U_hat/(Z'*M2*D/n);
       
end